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1  Introduction 


This  user’s  guide  presents  detailed  discussions  of  the  theoretical  aspects  of 
the  sigma  version  of  CH3D-WES  (Curvilinear  Hydrodynamics  in 
3-Dimensions  -  Waterways  Experiment  Station).  The  governing  equations  are 
presented  in  both  Cartesian  and  boundary-fitted  form  along  with  a  discussion 
of  boundary  conditions  and  solution  techniques.  Particular  attention  is  paid  to 
addressing  recent  improvements  to  the  representation  of  the  horizontal 
momentum  diffusion  terms  and  to  the  prediction  of  vertical  turbulent 
transport.  In  addition,  the  structure  of  the  computer  code  is  discussed  via  a 
description  of  the  order  and  function  of  each  subroutine.  Finally,  the  required 
input  to  operate  the  model  and  output  files  generated  by  the  model  are 
discussed. 
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2  Sigma  Stretched  CH3D- 
WES  Hydrodynamic  Model 


The  numerical  hydrodynamic  model  CH3D-WES  exists  in  both  a  Z-grid 
and  a  sigma  stretched  version  for  representation  of  the  vertical  dimension. 

The  Z-grid  version  was  developed  during  a  study  on  Chesapeake  Bay  and  is 
documented  in  Johnson,  et  al.  (1991b).  The  basic  sigma  stretched  model  was 
developed  by  Sheng  (1986)  for  WES  but  has  been  extensively  modified. 

These  modifications  have  consisted  of  implementing  different  basic  numerical 
formulations  of  the  governing  equations  as  well  as  substantial  recoding  of  the 
model  to  provide  more  efficient  computing.  In  particular,  two  recent 
modifications  presented  in  this  report  include  the  incorporation  of  a  compact 
form  of  the  horizontal  momentum  diffusion  terms,  and  a  two-equation  vertical 
(k-e)  turbulence  model.  As  its  name  implies,  CH3D-WES  makes  hydro- 
dynamic  computations  on  a  curvilinear  or  boundary-fitted  planform  grid. 
Physical  processes  impacting  circulation  and  vertical  mixing  that  are  modeled 
include  tides,  wind,  density  affects  (salinity  and  temperature),  freshwater 
inflows,  turbulence,  and  the  effect  of  the  earth’s  rotation. 

The  boundary-fitted  or  curvilinear  coordinate  feature  of  the  model  in  the 
horizontal  dimensions  provides  the  grid  resolution  enhancement  necessary  to 
adequately  represent  deep  navigation  channels  and  irregular  shoreline  configu¬ 
rations  of  the  flow  system.  The  curvilinear  grid  also  permits  adoption  of 
accurate  and  economical  grid  schematization  software.  The  solution  algorithm 
employs  both  an  external  mode,  consisting  of  vertically  averaged  equations 
which  provide  a  solution  for  the  free  surface  displacement  and  vertically 
averaged  velocities,  and  an  internal  mode.  The  deviation  of  the  horizontal 
components  of  the  full  3D  velocity  from  the  vertically-averaged  velocity 
components  are  computed  in  the  internal  mode  and  then  added  to  the  vertically 
averaged  components  to  yield  the  full  3D  horizontal  components.  In  addition, 
the  vertical  component  of  the  3D  velocity  field  and  the  3D  salinity  and 
temperature  fields  are  computed  in  the  internal  mode. 


Governing  Equations 

The  .governing  partial  differential  equations  are  based  on  the  following 
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assumptions:  a)  the  hydrostatic  pressure  distribution  adequately  describes  the 
vertical  distribution  of  fluid  pressure,  b)  the  Boussinesq  approximation  is 
appropriate,  c)  the  eddy  viscosity  approach  adequately  describes  turbulent 
mixing  in  the  flow. 

The  basic  equations  for  an  incompressible  fluid  in  a  right-handed  Cartesian 
coordinate  system  (x,y,z)  are  (Johnson  et  al.,  1991b): 


dx  dy  dz 


du  ^  du^  ^  duv  ^  duw  _  a,  _  ^  ^ 

*  ~dx  ^  ~dy  ^  Po  5^ 


du 
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dz 
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P  =  P  (T,S) 


(7) 


where 

iu,v,w)  =  velocities  in  (x,y,z)  directions 
t  =  time 

/  =  Coriolis  parameter  defined  as  2Q  sin  (j) 
where 

n  =  rotational  speed  of  the  earth 
(j)  =  latitude 
p  =  density 
p  =  pressure 

Ai^K,,  =  horizontal  turbulent  eddy  viscosity/diffusivity  coefficients 
A^,K^  =  vertical  turbulent  eddy  viscosity/diffusivity  coefficients 
g  =  gravitational  acceleration 
T  =  temperature 
S  =  salinity. 

Equation  4  implies  that  vertical  accelerations  are  negligible  and  thus  the 
pressure  is  hydrostatic.  Various  forms  of  the  equation  of  state  can  be 
specified  for  Equation  7.  In  the  present  model,  the  formulation  given  below 
is  used: 


p  =  P/(a  +  0.698P) 


(8) 


where 


p  =  density  in  grams  per  cubic  centimeter 
P  =  5890  +  38r  -  0.375r  +  35 
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a  =  1779.5  +  11.25r  -  0.07457’^  -  (3.8  +  0.017)5 


and  T  is  temperature  in  degrees  Celsius  and  S  is  salinity  in  parts  per  thousand 
(PPt)- 

Non-Dimensionalization  of  Equations 


The  dimensionless  forms  of  the  governing  equations  are  used  to  facilitate 
relative  magnitude  comparisons  of  the  various  terms  in  the  governing 
equations.  The  following  dimensionless  (star)  parameters: 


{u,  V,  w*)  =  («,  V,  >vxyz,)/c/, 
(x‘,  y‘,  z*)  =  (X,  y,  zXJZ^r 
{r\,r;)  =  {T\,T'^y)lpJZ,U, 

t*  =  tf 

r  =  gm^r  =  r/5. 

P*  =  (P-  Po)KPr  -  Po) 
r  =  (T-T,)I(T,-T,) 

A*  =  AJA„ 

K;  =  KJK,, 


K;  =  KJK,, 


where  (t\,  f^y)  =  wind  stress  in  (x,y)  direction,  p„  and  T„  are  the  expected 
minimum  density  and  temperature,  and  f  =  water  surface  elevation.  The  dim 
ensionless  parameters  resulting  from  these  definitions  are  the  Rossby  Number 
R„,  Froude  Number  F„  Densimetric  Froude  Number  and  the  horizontal 
and  vertical  Ekman  Numbers  E^,  E„  respectively.  These  dimensionless 
parameters  are  defined  as  follows: 


Ro  =  ujfx. 


Fr  =  UJigZrr^ 

Fr,  =  pJ'^UJ[gZXPr-P„)r  (10) 
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En  =  AJfX^ 

E.  =  AJfZ,^ 

All  parameters  with  r  subscripts  are  arbitrary  reference  scaling  quantities. 


Boundary-Fitted  Equations 

The  CH3D-WES  model  utilizes  a  boundary-fitted  or  generalized  curvilinear 
planform  grid  which  can  be  made  to  conform  to  flow  boundaries,  providing  a 
detailed  resolution  of  the  complex  horizontal  geometry  of  the  flow  system. 

This  necessitates  the  transformation  of  the  governing  equations  into  boundary- 
fitted  coordinates  (^,17).  If  only  the  (x,y)  coordinates  are  transformed,  a 
system  of  equations  similar  to  those  solved  by  Johnson  (1980)  for  vertically 
averaged  flow  fields  is  obtained.  However,  in  the  CH3D-WES  model  not 
only  are  the  (x,y)  coordinates  transformed  into  the  (I,??)  curvilinear  system  but 
the  velocity  also  is  transformed  such  that  its  components  are  contravariant 
(i.e.,  perpendicular  to  the  (^,17)  coordinate  lines).  This  is  accomplished  by 
employing  the  definitions  below  for  the  components  of  the  Cartesian  velocity 
(u,v)  in  terms  of  contravariant  components  u  and  v 


U  =  XM  +  X  V 

^  n 

V  =  y^u  +  y^v  (11) 


along  with  the  following  expressions  for  replacing  Cartesian  derivatives 


/,  =  \  WyX  -  (fy,)  ] 

^  J  ^  ^  T1 

fy.=  \  (^2) 

^  J  '  'I 


where /is  an  arbitrary  variable  and  J  is  the  Jacobian  of  the  coordinate 
transformation  defined  as 

J  =  xyy  -  X  y,  . 


(13) 


Additional  metric  coefficients  of  the  transformation  are: 
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✓t  2  2 


^22  =  xl  +  yl 

<^12  =  ^4^  ^  >^4^,  =  ^2V 


(14) 


In  addition  to  the  horizontal  grid  transformation,  the  vertical  dimension  is 
transformed  into  a  sigma-stretched  grid  (Figure  1)  by  : 

a  =  ^  ~  ^  (15) 

C  +  /i 

where  h  is  the  water  depth  from  the  datum  where  z  =  0. 


With  both  the  Cartesian  coordinates  and  velocity  components  transformed,  the 
following  non-dimensional  three  dimensional  governing  equations  are  solved 
on  the  sigma  stretched  grid. 
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dt  ^  J 


rr^CO 

+  H -  =  0 

dc 


(16) 


where 

H  is  the  total  water  depth,  i.e.,  (h+('). 


dHu 

~W 


22  9(^ 


Gn  9C^  ^ 

■  — -  — )  +  —  (G.m  +  G„v) 
dr]  y 
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^[^{Jy^Huu.JyHu^ 
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dx\ 


(Jy^Hu  V + Jy^Hv^] 


[-^  (Jx^Huu  +Jx^Hu 


+ 


d 


{Jx^Hu  V + Jx^HV^} 


-R  +  —  —  —)+ X  -  Horizontal  Diffusion 

f’  do  H  da  'da 


virt  °(^21  dp  <^12  9p 


Fr] 


{H\\:2L^-l^^)da 

J2  5^  jl 


■( 


^22  dH  _  <5,2  dH 
~jF-^  ~fF-^ 


(17) 


dHv 

dt 


-H{- 


G,  5C  - 


H  9^  H  dx\  J 


)--r((?l,M  +  <52i)^ 


r  d  ,r  rr—  r  rr— ^  9 
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-R  +  ^  —  (A  —)  +  Y- Horizontal  Diffusion 

°  da  H  da  '’da 


^^H 0,  5p  ^11  5p  • 


Fr: 


H  di,  F  dr\ 


^  n  dl  F  5n  Jc> 


ff 


(18) 


From  Equation  16,  the  vertical  velocity  in  the  sigma  coordinates,  w,  is 
computed,  i.e.. 


fi>. 


top 


%or 


IT  ^  J 


■W  (■'^"1 

d^ 


+ 


dr] 


1 


(19) 


where  cj,„p  is  the  vertical  velocity  at  the  top  of  a  cell  and  is  the  velocity  at 
the  bottom  of  the  cell.  The  computations  proceed  from  the  bottom  of  the 
water  column,  where  =  o,  to  the  surface. 

The  vertical  velocity  of  a  water  particle,  w,  at  some  a  location  can  then  be 
computed  from 


w  =  H(i)  + 


1+a  S(^ 

"p" 


(20) 


where 

p  =  gzjl/xf) 


The  horizontal  diffusion  terms  in  Equations  17  and  18  are  rather  lengthy  and 
as  a  result  are  presented  separately  in  Appendix  A.  The  equations  shown  in 
Appendix  A  are  a  compact  form  of  those  derived  by  Johnson  et  al.  1991b.  A 
description  of  the  implementation  of  the  horizontal  diffusion  terms  is  provided 
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in  Chapman  (1993).  The  transport  equations  for  salt  and  temperature  are 
written; 


dHS  _  Q  dS^  _  p  dHaS  _  Ec  ( dJHuS  ^  dJHvS^ 

dt  Pr^H  dc  "  8a  ‘’da  J  9^  ^'*1 


8 

^22 

K.H  -Ji 

8S 

9 
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'  J 
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<^12 
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8S 
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*  J 

9p 

J  9r| 

dHT  ^  E,  8  dT.  _  ^  8H&T  _  .  dJHuT  ^  dJHvT. 
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Gjj 

KH 

8T 
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_94 

94  J 

*  Jdr\ 

8 

G,2 

K.H  — li 

8T 

9 

K,H^  E. 

9ti 

94  J 

9p 

*  J  9q 

Pr,  and  Pr^  are  vertical  and  horizontal  Prandtl  numbers  which  are  ratios  of 
eddy  viscosity,  or  A^,  to  eddy  diffiisivity,  and  K^. 

Similarly,  the  transformed  external  mode  equations  are  written: 


dt 


J  9^ 


(23) 


8t 


^  9C 
7  9r| 
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^^{l-{Jy,UU^JyUV)*^{Jy,UV*JyVV)] 

jm  dC  ’  dr\  ^ 

-  ^  [1.  {JxUU-^Jxm  +  A  (Jx.UV^Jx  F  jO] 

J^H  dV  '  ^  5n  ' 

+  +  X  -  Horizontal  Diffusion 

-  ^^1(—  A  -  —  -A) 

2Frj  -P  P  dr] 


and 


A  H( 

5r  P 


^[^(Jy^UU*JyUV)*±{Jy^UV*JyVr)] 
PH  d%  ^  dr\  ^ 


+  ^[^iJxMU-rJx  UV)  +  —  (Jx.UV+JxVV)] 
PH  dC  ^  '  Sri  ^ 


+  -  x^^  +  y  -  Horizontal  Diffusion 


2Fr]  P  di,  P  9ti 


Boundary  Conditions 


The  boundary  conditions  at  the  free  surface  are 


A 


du  dv 
dz  ’  dz 


,  x_yp  .{civl.c  (fj) 


^  ~  —  K  (T  - 

dz  E„ 


(25) 


(26) 
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whereas  the  boundary  conditions  at  the  bottom  are 


'  du  dv' 
dz  ’  dz 


U  _,\l/2  _  _ 

=  ^ZC,(„f*vf)  V,) 

vr 

dT  dS  ^  ^ 

■&  ’  ■& 


(27) 


where 

Tj  =  wind  shear  stress 
C  =  surface  drag  coefficient 
W  =  wind  speed 

K  =  surface  heat  exchange  coefficient 

Te  =  equilibrium  temperature 

Tb  =  bottom  shear  stress 

Q  =  bottom  friction  coefficient 

u,  ,  Vi  =  near  bottom  horizontal  velocity  components. 

Defining  Z;  as  one  half  of  the  bottom  layer  thickness  and  assuming  a  log 
velocity  profile,  Q  is  given  by 

C,  =  k^[\n{zjz,)r  (28) 

where  k  is  the  von  Karman  constant  (0.4)  and  Zq  is  a  bottom  roughness  height. 
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Manning’s  formulation  is  employed  for  the  bottom  friction  in  the  external 
mode  equations  if  the  model  is  used  only  to  compute  vertically  averaged  flow 
fields.  The  surface  drag  coefficient  is  computed  according  to  Garratt  (1977) 
as  follows: 

C  =  (0.75  +  0.067  JV)  X  10'^ 


with  the  maximum  allowable  value  being  0.003.  The  surface  heat  exchange 
coefficient,  K,  and  the  equilibrium  temperature,  T„  are  computed  from 
meteorological  data  (wind  speed,  cloud  cover,  wet  and  dry  bulb  air  tempera¬ 
tures,  and  relative  humidity)  as  discussed  by  Edinger,  Brady,  and  Geyer 
(1974).  The  wind  speed,  W,  must  be  in  meters/second. 

Freshwater  inflow  and  water  temperature  are  prescribed  along  the  shoreline 
where  river  inflow  occurs,  however,  the  salinity  at  the  river  boundary  is 
specified  according  to  a  zero  spatial  gradient  assumption  (computed  from  the 
previous  time  step).  At  an  ocean  boundary,  the  water-surface  elevation  is  pre¬ 
scribed  along  with  time-varying  vertical  distributions  of  salinity  and  tempera¬ 
ture.  Specified  values  of  salinity  and  temperature  are  employed  during  flood 
flow,  whereas,  during  ebb,  interior  values  are  advected  out  of  the  grid.  The 
normal  component  of  the  velocity  and  the  eddy  viscosity  and  diffusivity  are 
set  to  zero  along  solid  boundaries. 


Initial  Conditions 


When  initiating  a  run  of  CH3D-WES,  the  values  of  f,  m,  v,  w,  U  and  V 
are  set  to  zero.  Values  of  salinity  and  temperature  are  read  from  input  files. 
These  initial  data  are  generated  from  prototype  measurements  at  a  limited 
number  of  locations.  Once  the  values  in  individual  cells  are  determined  by 
interpolating  from  the  field  data,  the  resulting  3D  field  is  smoothed. 
Generally,  the  salinity  and  temperature  fields  are  held  constant  for  the  first 
few  days  of  a  simulation. 


Computational  Grid 


A  staggered  grid  is  used  in  both  the  horizontal  and  vertical  directions  of 
the  computational  domain  (Figure  2).  In  the  horizontal  direction,  a  unit  cell 
consists  of  a  f-point  in  the  center  (^y),  a  U-point  to  its  left  (t/y),  and  a  V- 
point  to  its  bottom  (V^y).  In  the  vertical  direction,  the  vertical  velocities  are 
computed  at  the  "full"  grid  points.  Horizontal  velocities,  temperature, 
salinity,  density  and  turbulence  quantities  are  computed  at  the  "half  grid 
points  (half  grid  spacing  below  the  full  points). 

Two  arrays,  each  of  dimension  (IMAX,  JMAX),  are  used  to  index  the  grid 
cells.  The  array  NS  indicates  the  condition  of  the  left  and  right  cell 
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Figure  2.  Staggered  grid 


boundaries,  while  the  array  MS  denotes  the  condition  of  the  top  and  bottom 
cell  boundaries  (Figure  3). 


Numerical  Solution  Algorithm 

Finite  differences  are  used  to  replace  derivatives  in  the  governing  equa¬ 
tions,  resulting  in  a  system  of  linear  algebraic  equations  to  be  solved  in  both 
the  external  and  internal  modes.  The  external  mode  solution  consists  of  the 
surface  displacement  and  vertically  integrated  contravariant  unit  flows  U  and 
V.  All  of  the  terms  in  the  transformed  vertically-averaged  continuity  equation 
are  treated  implicitly,  whereas,  only  the  water-surface  slope  terms  in  the  trans¬ 
formed  vertically-averaged  momentum  equations  are  treated  implicitly.  If  the 
external  mode  is  used  as  purely  a  vertically-averaged  model,  the  bottom  fric¬ 
tion  is  also  treated  implicitly.  Those  terms  treated  implicitly  are  weighted 
between  the  new  and  old  time-steps.  The  resulting  finite  difference  equations 
are  then  factored  such  that  a  |-sweep  followed  by  an  rj-sweep  of  the  horizontal 
grid  yields  the  solution  at  the  new  time-step.  Writing  Equations  23-25  as 
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PA? 

AriJ 


^V)'  -  Ur); 


and 


QAtHG^^  I  ^ 


(33) 


(34) 


The  Tj-sweep  then  provides  the  updated  f  and  V  at  the  n  + 1  time  level. 


r|  -sweep 


p9A? 

Ar|J 


C,;  -  (i-e) 


pA? 

A^ 


(35) 


and 


^+l 


6A fUG^^  /  «+i  «+i\ 

-  Q- ) 


(1-0) 


A^ 


C”)  +  A?  N 


(36) 


A  typical  value  of  6  of  0.55  is  employed.  M  and  N  represent  all  terms  in  the 
equations  evaluated  at  the  previous  time  step. 
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The  internal  mode  consists  of  computations  for  the  three  velocity  compo¬ 
nents  u,  V,  and  w,  salinity,  and  temperature.  Defining  the  horizontal 

components  of  the  3D  velocity  as  u  =  IL  +  u'  and  v  =  ^  +  v',  the 

differential  equations  for  the  {u',  v' )  components  are  obtained  by  subtracting 
the  vertically  averaged  momentum  equations  (23-25)  from  the  3D  momentum 
equations  (17-18).  This  removes  the  water  surface  slope  terms  from  the  equa¬ 
tions  of  motion  which  removes  the  restrictive  free-surface  gravity  wave  speed 
from  the  internal  mode  stability  criteria. 

It  is  important  to  ensure  that  the  vertical  integration  of  the  («',  v')  is  zero. 
This  is  accomplished  by  evaluating  the  nonlinear  inertia  and  turbulent  diffu¬ 
sion  terms  in  the  vertically-averaged  momentum  equations  by  summing_  the 
corresponding  terms  in  the  3D  equations  at  all  vertical  layers.  Once  (m  ,  v  ) 
are  determined,  they  are  slightly  adjusted  to  absolutely  ensure  that  their  verti¬ 
cal  sum  is  zero  and  then  are  added  to  the  vertically  averaged  velocities  to 
yield  the  horizontal  components  of  the  full  3D  velocity. 

The  only  terms  treated  implicitly  are  the  vertical  diffusion  terms  in  all 
equations  and  the  bottom  friction  in  the  momentum  equations.  Roache’s 
(1972)  second  upwind  differencing  is  used  to  represent  the  convective  terms  in 
the  momentum  equations,  whereas,  a  spatially  and  temporally  third-order 
scheme  developed  by  Leonard  (1979)  called  QUICKEST  is  used  to  represent 
the  advective  terms  in  Equations  21  and  22  for  salinity  and  temperature,  re¬ 
spectively.  For  example,  if  the  velocity  on  the  right  face  of  a  computational 
cell  is  positive  then  the  QUICKEST  value  of  the  salinity  computed  for  the  flux 
through  the  face  is 


U.  , At 


(37) 


The  more  interested  reader  is  referred  to  the  paper  by  Leonard  (1979). 


Two  Equation  k  -  e  Turbulence  Closure 

A  vertical  k-e  turbulent  eddy  viscosity  model  which  includes  the  effects  of 
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wind  shear,  bottom  shear,  velocity  gradient  turbulence  production,  dissipation, 
diffusion  and  stratification  has  been  implemented.  The  basic  idea  behind  the 
(k-e)  turbulence  model  (Rodi,  1980;  ASCE,  1988)  is  that  the  vertical  eddy 
viscosity  coefficient  can  be  related  to  the  turbulence  energy  per  unit  mass,  t, 
its  rate  of  dissipation,  e;  and  an  empirical  coefficient  (c,  =  0.09),  i.e.: 


A 


z 


(38) 


The  transport  equations  for  the  turbulence  quantities  are  written: 


d{Hk)  _  1  L 
dt  H  da  ^  da 


(P  -  8  +  G)H 


(39) 


and 


djHs) 

dt 


1  d  ds 

H  da  a  da 

£ 


=  (c,iP  -c,—)H 
k  "  k 


(40) 


in  which  =  1.3,  C;  =  1.44,  and  C2  =  1.92.  The  source  and  sink  terms  on 
the  right  hand  side  of  equations  39  and  40  represent  mechanical  production  of 
turbulence,  P^,  due  to  vertical  velocity  gradients  and  buoyancy  production  or 
dissipation,  G.  The  functional  forms  of  these  mechanisms  are  as  follows: 


II 

I 

I 

du 

'  +  2G,2 

du  dv 

'  dv' 

2 

^  IP 

11 

do 

da  da 
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da 

and 


G  =  S 
HPr^  p  da 


(42) 


where,  as  previously  noted,  is  the  turbulent  Prandtl  Number.  Surface  and 
bottom  boundary  conditions  for  the  turbulence  quantities  are  specified  as 
follows: 
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(43) 


and 


U\ 

WKg 


(44) 


where  k  is  the  von  Karman  coefficient.  The  friction  velocity,  C/.,  at  the  sur¬ 
face  boundary  is  defined  as  the  square  root  of  the  resultant  wind  shear  stress, 
TSR,  where: 


TSR  -  ^ 


(45) 


in  which  and  Ty  are  the  components  of  the  wind  stress.  The  bottom  friction 
velocity  is  computed  in  an  identical  way  with  the  wind  vectors  replaced  by  the 
contravariant  velocity  components.  The  suppression  of  the  vertical  diffusivity 
by  stratification  is  accomplished  by  modifying  the  computed  value  as  follows: 


where  /?,  is  the  Richardson’s  Number  (Bloss  et  al.  1988). 


g  dp 

pH  da 

1 

2 

1 

d\u^  + 

Ip 

da 

(47) 


A  complete  description  of  the  implementation  of  the  k-c  model  is  presented  in 
Chapman  (1994). 
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3  Structure  of  Sigma 
CH3D-WES 


The  CH3D-WES  model  has  a  main  program  as  well  as  several  subroutines. 
Subroutines  governing  model  setup  are  called  from  the  main  program  while 
subroutines  governing  computations  are  called  from  subroutine  CH3DM2. 
Each  of  these  is  listed  below  with  a  description  of  its  function.  Entry  points 
in  subroutines  are  also  noted.  Two  INCLUDE  files,  application. inc  and 
ch3d.inc,  are  needed.  They  are  used  to  set  up  parameters,  dimensions  of 
various  arrays,  and  COMMON  blocks.  During  model  compilation,  these  files 
are  inserted  wherever  the  INCLUDE  statements  in  the  source  code  call  for  the 
files.  Several  input  data  files  are  required.  These  are  listed  in  Appendix  C. 


CH3D  The  main  program. 

CH3DIR  Reads  data  from  main  input  file,  FILE  4  (see  Appendix  B), 

which  controls  computations,  input,  and  output.  Various  con¬ 
stants  are  computed,  and  the  vertical  (a-)  layer  thicknesses  are 
set. 


CH3DTR  Reads  (x,y)  coordinates  (ft)  and  depths  (ft)  at  the  cell  corners 
of  the  boundary-fitted  grid  from  FILE  15  (ITRAN=2).  The 
coordinates  are  then  multiplied  by  the  scale  factor,  XMAP, 
and  divided  by  XREF  to  make  them  nondimensional. 
Subroutine  BJINTR  is  called  to  provide  the  coordinate  deriva¬ 
tives  needed  to  compute  the  metrics  of  the  transformation. 

BJINTR  Computes  various  coordinate  derivatives  and  sets  the  water 

depths  HU(I,J)  and  HV(I,J)  on  the  faces  of  each  computational 
cell. 


CH3DIH  Prints  water  depths,  if  requested  by  input  data.  Also,  the 

water  depths  are  made  nondimensional  by  dividing  by  ZREF. 


CH3DND  Normalizes  several  variables  and  parameters,  such  as  the 
Ekman  number,  Rossby  number,  time-step,  etc. 
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CH3DII  Sets  up  the  arrays  of  boundary  flags  that  indicate  the  nature  of 

computational  cell  boundaries.  In  addition,  arrays  controlling 
the  computation  of  the  convective  terms  in  the  momentum 
equations  and  the  water  surface  cross-derivative  terms  are  set 
up.  One-dimensional  channel  cells  are  identified. 

CH3DIF  Initializes  various  variables  for  a  cold  start  run  and  opens  time 
series  output  files  for  elevation,  velocities,  salinity,  etc.  as 
well  as  print  and  snapshot  files.  The  hot  start  capability  is  not 
operational. 

CH3DIV  The  arrays  created  in  CH3DII  concerning  water  surface  cross¬ 
derivatives  contain  logical  values.  Those  arrays  are  used  in 
this  subroutine  to  create  arrays  containing  numerical  values. 
These  arrays,  i.e.  AFV1(I,J),  etc.  are  used  to  control  compu¬ 
tation  of  not  only  water  surface  cross-derivatives  but  other 
variables  as  well. 

CH3DWS  Controls  the  reading  of  either  wind  speed  or  wind  stress.  If 
the  wind  speed  is  read,  the  stress  is  computed  from  Garratt’s 
equation.  ENTRY  CH3DWT  controls  the  time-varying  reads 
and  computations. 

The  subroutines  above  are  called  from  CH3D  in  the  sequence  given. 

Before  calling  CH3DM2,  which  controls  the  computations,  the  initial  salinity 

field  is  read  from  FILE  74.  The  initial  temperature  field  is  read  from 

FILE  17  and  made  dimensionless. 

CH3DM2  Final  subroutine  called  fi:om  CH3D.  All  subroutines  control¬ 
ling  the  actual  3D  computations  are  called  from  this  subroutine 
in  the  order  they  appear  below. 

CH3DDP  Computes  the  total  water  depths  from  the  latest  water  surface 
elevation  field.  ENTRY  CH3DDM  sets  total  water  depths  at 
the  intermediate  time  level  M  and  ENTRY  CH3DDN  sets  total 
water  depths  at  time  level  N. 

CH3DTK  Reads  equilibrium  temperatures  and  surface  heat  exchange 

coefficients  from  FILE  19  and  then  casts  them  into  nondimen- 
sional  form.  ENTRY  CH3DTB  controls  the  time-varying  read 
and  interpolation. 

CH3DRI  Reads  river  inflows  from  FILE  13.  ENTRY  CH3DRV  con¬ 
trols  the  time-varying  read  and  interpolation. 

CH3DTI  Reads  and  initializes  tidal  boundary  conditions  from  FILE  16. 
ENTRY  CH3DTD  updates  boundary  values. 
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CH3DSAI 

CH3DTEI 

CH3DDE 

CH3DKE 

CH3DWT 

CH3DTB 

CH3DRV 

CH3DTEV 

CH3DTD 

CH3DSAV 

CH3DDN 

CH2DXY 

CH3DDP 

CH3DXYZ 

CH3DDI 

CH3DSA 


Reads  salinity  and  temperatures  at  tidal  boundaries  from  FILE 
76.  ENTRY  CH3DSAV  controls  the  time-varying  reads, 
interpolation,  and  conversion  to  nondimensional  form. 

Reads  temperatures  at  river  inflow  boundaries  from  FILE  78. 
ENTRY  CH3DTEV  controls  the  reading  of  time-varying 
temperatures  and  interpolation. 

Computes  the  water  densities  using  Equation  (8).  The  baro- 
clinic  terms  in  the  momentum  equations  are  then  evaluated. 

Computes  the  eddy  viscosity  and  eddy  diffiisivity  coefficients. 

ENTRY  in  CH3DWS  for  reading  time-varying  wind  data. 

ENTRY  in  CH3DTK  for  reading  time-varying  equilibrium 
temperature  and  heat-exchange  coefficient. 

ENTRY  in  CH3DRI  for  reading  time- varying  river  flows. 

ENTRY  in  CH3DTEI  for  reading  time-varying  temperatures  at 
river  inflow  boundaries. 

ENTRY  in  CH3DTI  for  reading  time-varying  tide  data. 

ENTRY  in  CH3DSAI  for  reading  time-varying  salinities  and 
temperatures  at  tidal  boundaries. 

ENTRY  in  CH3DDP  for  assigning  total  water  depths  at  time 
level  N. 

Computes  the  vertically  averaged  flow  field  from  the  vertically 
averaged  equations  of  motion. 

Using  the  water  surface  field  computed  in  CH2DXY,  com¬ 
putes  the  total  water  depths  at  time  level  N  -H 1 . 

Computes  the  3D  velocity  field.  Mass  conservation  is  ensured 
by  forcing  the  vertical  sum  of  the  horizontal  components  of 
the  3D  velocity  to  match  the  vertically  integrated  values 
computed  in  CH2DXY. 

Computes  the  convective  and  diffusion  terms  in  the  momentum 
equations  using  the  most  recent  computation  results  from 
CH3DDP  and  CH3DXYZ.  These  terms  are  then  employed  at 
the  next  time  step  in  CH2DXY  and  CH3DXYZ. 

Computes  the  salinity  field. 
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CH3DTE 


Computes  the  temperature  field. 

CH3DBL  Checks  the  water  surface  elevations  for  the  program  "blowing 
up". 

CH3DOT  Controls  the  output  printed  and/or  written  to  files  for  plotting. 

Output  is  in  terms  of  physical  dimensional  variables.  Subrou¬ 
tine  CH3DC1  is  called  with  ENTRIES  CH3DC2,  CH3DC3, 
CH3DC4,  CH3DC5,  CH3DC6,  CH3DC7,  CH3DC8, 
CH3DC9,  CH3DCA,  CH3DCC,  CH3DCD,  and  CH3DCE. 
Each  is  described  below. 

CH3DC1  Provides  dimensional  water  surface  elevations. 

CH3DC2  Provides  dimensional  physical  vertically-averaged  velocity  in 
x-direction. 

CH3DC3  Provides  dimensional  physical  vertically-averaged  velocity  in 
y-direction. 

CH3DC4  Provides  dimensional  physical  horizontal  velocity  component 
in  x-direction. 

CH3DC5  Provides  dimensional  physical  horizontal  velocity  component 
in  y-direction. 

CH3DC6  Provides  dimensional  physical  vertical  component  of  3D 
velocity. 

CH3DC7  Provides  salinity. 

CH3DC8  Provides  dimensional  temperature. 

CH3DC9  Provides  dimensional  physical  magnitude  and  direction  of 

horizontal  velocity. 

CH3DCA  Provides  dimensional  physical  horizontal  components  of  3D 
velocity  at  the  centers  of  cells. 

CH3DCC  Provides  dimensional  water  density. 

CH3DCD  Provides  dimensional  vertical  eddy  viscosity. 

CH3DCE  Provides  dimensional  vertical  eddy  diffiisivity. 

In  subroutine  CH3DOT,  the  following  files  are  created  for  use  in 

generating  time  series  plots,  vector  plots,  or  contour  plots. 
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FILE  21  For  time  series  plots  of  dimensional  water  surface  elevation  at 
specified  horizontal  locations. 

FILE  22  For  time  series  plots  of  dimensional,  Cartesian  horizontal 

velocities  (x  and  y  directions)  at  cell  centers  at  specified  hori¬ 
zontal  locations. 

FILE  23  Geometry  of  study  area  (needed  for  plotting  snapshots  or  con¬ 
tours). 

FILE  24  For  velocity  vector  plots  and  contour  plots  of  surface 
elevation,  salinity,  temperature,  etc. 

FILE  25  For  time  series  plots  of  discharges  at  specified  horizontal 
ranges. 

FILE  31  For  time  series  plots  of  salinity  at  specified  horizontal  loca¬ 
tions  in  all  layers. 

FILE  34  For  time  series  plots  of  temperature  at  specified  horizontal 
locations  in  all  layers. 

FILE  35  For  time  series  plots  of  vertical  eddy  viscosity  at  specified 
horizontal  locations  in  all  layers. 

FILE  36  For  time  series  plots  of  vertical  eddy  diffusivity  at  specified 
horizontal  locations  in  all  layers. 

FILE  37  For  time  series  plots  of  density  at  specified  horizontal 
locations  in  all  layers. 


As  previously  indicated,  there  are  two  INCLUDE  files,  application. inc  and 
ch3d.inc,  needed  for  running  the  model.  Of  these,  application. inc  is  used  to 
set  up  the  model  parameters  for  running  the  particular  application.  The  fol¬ 
lowing  parameters  are  set.  They  are  used  to  dimension  arrays  in  COMMON 
blocks  in  ch3d.inc  and  other  arrays  in  the  model.  Of  these,  ICELLS, 
JCELLS,  UMAX,  and  KM  have  to  be  set  exactly.  The  others  can  be  greater 
than  or  equal  to  what  is  needed.  The  ch3d.inc  file  does  not  have  to  be 
changed  from  application  to  application,  but  remains  the  same. 


ICELLS 

JCELLS 

UMAX 

KM 

NSTATS 


:  Number  of  grid  cells  in  the  ^-direction 
:  Number  of  grid  cells  in  the  17-direction 
;  The  greater  of  ICELLS  and  JCELLS,  plus  1 
:  Number  of  u-layers  in  the  vertical 

:  Maximum  number  of  gage  stations  where  information  will  be 
saved 


NTIDES  =  11  :  Not  used 

NRIVRS  :  Number  of  river  boundaries  used 
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NCNST  =  37 


NBNDS 

NBARRS 

NPRWIN 

NSNAPS 

NRANGS 

NTIDFN 

NTIDBN 

NTIDPT 

NROWS 

NCOLS 

KROWS 

NX8PTS 

NY8PTS 

SPVAL 


:  Maximum  number  of  tidal  constituents  used  (set  to  37)  - 
used  only  if  tidal  signals  were  generated  using  constituents 
not  operational. 

:  Number  of  open  water  boundaries  used 
:  Number  of  interior  thin-wall  barriers 
:  Number  of  print  windows  for  printing  model  results 
:  Number  of  snapshot  windows  where  information  is  saved 
:  Number  of  ranges  where  discharge  information  is  saved 
:  Number  of  tide  ftinctions  used 
:  Number  of  tidal  boundaries  used 
;  Maximum  number  of  values  in  the  input  tide  functions 
;  Maximum  number  of  computational  chains  used  in 
^-direction 

:  Maximum  number  of  computational  chains  used  in 
ij -direction 

;  Larger  of  NROWS  and  NCOLS,  plus  1 
:  Number  of  one-cell  wide  channel  cells  in  ^-direction 
:  Number  of  one-cell  wide  channel  cells  in  rj-direction 
:  A  small  value  to  which  the  vertical  eddy  coefficients,  etc. 
are  set  as  a  default 
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4  Summary 


The  purpose  of  this  report  is  to  describe  the  main  features  of  the  sigma 
version  of  the  CH3D-WES  hydrodynamic  model.  In  Chapter  2,  the  basic 
governing  equations  are  given  followed  by  the  boundary  and  initial  conditions 
employed.  The  report  outlines  the  structure  of  the  computer  model  in 
Chapter  3,  listing  the  names  of  the  various  subroutines,  their  functions,  and 
the  calling  sequence.  This  should  help  users  who  are  interested  in  following 
the  logic  of  the  model.  Also  listed  are  various  output  files  created  by  the 
model  and  their  contents. 
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Appendix  A 
Transformed  Horizontal 
Momentum  Diffusion  Terms 


X  -  Horizontal  Diffusion 

=  Il  [  (XHu)  ;  +  (X^Hv)  ;] 

[  J  Jf 

+  ^  [  (XfHu)  +  (X,Hv)  ,] 

J  Jtj 

J  Jtj 

-  [:^  [  (XjHu)  ,  +  (X,Hv)  ,]  ' 

-  [  (X.Hu)  5  +  (X,Hv)  j] 

^  Jtj 

+  ^  [(y^Hu),+ (y,Hv),]  ^ 

+  ^  [  ( Y^Hu)  J  +  ( y,Hv)  j] 
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Y  -  Horizontal  Diffusion 

J  ^  U 

^  I  h 

+  ^  [  (  Y^HV)  J  +  ( Y^Hu)  j] 

u  ^ 

-  4i  ?] 

lJ  I  l7  ^ 

[  (i;Hv)f  +  (yjHu)f] 

I  ^  J^ 

+  ^4  [  (^„Hv)  ,  +  iX.Hu)  ,] 

I  iJ  ^  J^ 

Replacing  Hu  and  Hv  with  0  and  V,  respectively,  the  same  expressions  apply 
in  the  external  mode  equations. 
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Appendix  B 

List  of  Input  Data  in  Fiie  4 


DUMMY 

TITLE  Run  description  (Format  A80) 


DUMMY 

ITl,  IT2,  DT,  ISTART,  ITEST,  ITSALT  (2I8,F8.0,4I8) 

ITl  ;  Starting  time  step  (  always  set  =  1) 

IT2  ;  Ending  time  step 

DT  ;  Computational  time  step  in  sec 

ISTART  =  0  ;  Cold  start 

>0  ;  Hot  start  (not  operational) 

ITEST  =0  ;  No  diagnostic  output 

>  0  ;  Diagnostic  output 

ITSALT  ;  Number  of  time  steps  after  which  salinity  and 
temperature  computations  are  initiated 


DUMMY 

WPRCRD  (918, A8) 

WPRCRD  ;  Number  of  print  control  cards  which  follow 
DUMMY 

WXCELl,  WXCEL2,  WYCELl,  WYCEL2,  WZCELl,  WZCEL2,  WPRINT, 
WPRSTR,  WPREND,  WPRVAR  (918, A8) 


If  WPRCRD  >  0,  WPRCRD  cards  have  to  be  furnished  below. 


WXCELl 

WXCEL2 

WYCELl 

WYCEL2 

WZCELl 

WZCEL2 

WPRINT 

WPRSTR 

WPREND 

WPRVAR 


;  Starting  ^-cell  index 

;  Ending  ^-cell  index 

;  Starting  ij-cell  index 

;  Ending  ij-cell  index 

;  Starting  sigma  layer  index 

;  Ending  sigma  layer  index 

;  Printing  interval  (number  of  time  steps) 

;  Time  step  when  printing  starts 
;  Time  step  when  printing  ends 
;  Character  string  indicating  variables  printed 
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Note  :  The  following  characters  are  used  in  WPRVAR  for  designating 
different  variables. 

E  :  Surface  elevation  (cm) 

X  :  X-direction  unit  flow  rate  (cm^/sec) 

Y  :  Y-direction  unit  flow  rate  (cmVsec) 

U  :  X-direction  velocity  (cm/sec) 

V  :  Y-direction  velocity  (cm/sec) 

W  :  Z-direction  velocity  (cm/sec) 

S  :  Salinity  (ppt) 

T  :  Temperature  (deg  C) 

A  :  Average  velocity  magnitude  (cm/sec)  and  direction 
(measured  clockwise  from  the  true  North,  deg) 


DUMMY 
SNPCRD  (918, A8) 

SNPCRD  ;  Number  of  snapshot  control  cards  to  follow 
DUMMY 

SXCELl,  SXCEL2,  SYCELl,  SYCEL2,  SZCELl,  SZCEL2,  SNPINT, 
SNPSTR,  SNPEND,  SNPVAR  (9I8,A8) 

If  SNPCRD  >  0,  SNPCRD  cards  have  to  be  furnished  below. 


SXCELl 

SXCEL2 

SYCELl 

SYCEL2 

SZCELl 

SZCEL2 

SNPINT 

SNPSTR 

SNPEND 

SNPVAR 


Starting  ^-cell  index 

Ending  ^-cell  index 

Starting  r;-cell  index 

Ending  r?-cell  index 

Starting  sigma  layer  index 

Ending  sigma  layer  index 

Snapshot  interval  (number  of  time  steps) 

Time  step  when  snapshots  start 

Time  step  when  snapshots  end 

Character  string  indicating  snapshot 

variables  (same  notation  is  used  as  in 

WPRVAR) 


DUMMY 
NRANG  (9I8,A8) 

NRANG  ;  Number  of  ranges  for  computing  discharges 
DUMMY 

RANGDR,  RPOSl,  RPOS2,  RPOS3,  RRNAME  (7X,A1,3I8,A45) 
If  NRANG  >  0,  NRANG  cards  have  to  be  furnished  below. 


RANGDR 

RPOSl 

RPOS2 

RPOS3 

RRNAME 


Range  direction  (X  for  ^  and  Y  for  17) 
^  {rj)  cell  index  of  range  line 
Starting  t/  (^)  cell  index  for  range 
Ending  77  (0  cell  index  for  range 
Range  description  (name) 


DUMMY 
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IGI,  IGH,  IGT,  IGS,  IGU,  IGW,  IGC,  IGQ,  IGP  (1018)  :  Printout  flags, 
value  of  1  turns  printing  on  and  0  turns  it  off. 

IGI  ;  Print  arrays  such  as  NS,  MS,  NR,  MR,  etc. 


IGH 

IGT  =  0 
IGS  =  0 
IGU  =  0 
IGW  =  0 
IGC 


Print  all  depth  arrays 
Not  used 

Print  restart  arrays 
Not  used 
Not  used 

Print  grid  coordinates  and  depths 


IGQ  =  0 
IGP  ;  Save 


;  Not  used 

grid  information  in  FILE  23  for  plotting 
snapshots 


A 


DUMMY 

XREF,  ZREF,  UREF,  COR,  GR,  ROO,  ROR,  TO,  TR  (10F8.0) 


XREF 

;  Reference  horizontal  grid  distance 

(Maximum  horizontal  dimension  divided  by 
number  of  cells  in  that  direction,  cm) 

ZREF 

;  Reference  depth  (average  depth  in  cm) 

UREF 

;  Reference  horizontal  velocity 

(average  velocity  in  cm/sec) 

COR 

;  Coriolis  parameter 

GR 

;  Gravitational  acceleration  (cm/sec^) 

ROO 

;  Minimum  density  expected  (gm/cc) 

ROR 

;  Reference  density  (maximum  expected)  (gm/cc) 

TO 

;  Minimum  temperature  (Celsius) 

TR 

;  Reference  temperature  (maximum  expected) 
(Celsius) 

DUMMY 

THETA  (10F8.0) 

THETA  ;  Time  level  weighting  factor  in  computations 
(0.5  <  THETA  <  1.0) 


DUMMY 

ITEMP,  ISALT,  ICC,  IFI,  IFA,  IFB,  IFC,  IFD  (1018) 

ITEMP  =0  ;  No  computation  of  temperature 

=  -1  ;  Compute  temperature  (use  daily  equilibrium 
temperature  as  river  boundary  temperature) 

=  -2  ;  Compute  temperature  (use  time-varying 

temperature  as  river  boundary  temperature) 
ISALT  =0  ;  No  computation  of  salinity 

=  -2  ;  Compute  salinity,  setting  salinity  and 
temperature  at  tidal  boundaries 


ICC 

=  0 

;  Not  used 

IFI 

=  1 

Compute  nonlinear  (inertia)  terms 

= 

0 

No  computation  of  nonlinear  terms 

IFA 

=  0 

;  Not  used 

IFB 

=  0 

;  Not  used 
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IFC  =  0  ;  Not  used 

IFD  =  1  ;  Compute  horizontal  diffusion  terms 

=  0  ;  No  computation  of  horizontal  diffusion  terms 

DUMMY 

TWE,  TWH,  FKB  (3F8.0) 

TWE  ;  Temperature  in  the  epilimnion  (for  computing 
initial  conditions) 

TWH  ;  Temperature  in  the  hypolimnion  (for  computing 
initial  conditions) 

FKB  ;  Vertical  grid  index  of  the  initial 

thermocline  location  (for  computing  initial 
conditions) 

Note:  The  initial  conditions  computed  using  TWE,  TWH,  and  FKB  are 
overridden  by  FILE  17. 


DUMMY 

lEXP,  lAV,  AVR,  AVI,  AV2,  AVM,  AVMl,  AHR  (2I8,8F8.0) 


lEXP 

Vertical  eddy  coefficient  flag 

lEXP  =  0 

Constant  eddy  coefficient. 

=  1 

k  -  6  turbulence  closure 

lAV  =  0 

Not  used 

AVR 

Reference  vertical  eddy  viscosity  (cmVsec) 

AVI 

Not  used 

AV2 

Not  used 

AVMl  ; 

Minimum  allowable  vertical  eddy  diffusivity 
(cmVsec) 

AHR  ; 

;  Reference  horizontal  eddy  viscosity  or 

diffiisivity  (cmVsec) 


DUMMY 

GAMAX,  GBMAX  (2F8.0) 

GAMAX  ;  Maximum  value  of  eddy  viscosity  (cm^/sec) 
GBMAX  ;  Maximum  value  of  eddy  diffusivity  (cmVsec) 


DUMMY 
IWIND,  TAUX, 
IWIND  =  0 
=  1  ; 
=  2  ; 
=  3  ; 
TAUX  ; 


TAUY 


TAUY  (I8,5F8.0) 

;  Steady  and  uniform  wind  stress 
Steady  and  uniform  wind  speed 
Time  variable  and  uniform  wind  stress 
Time  variable  and  uniform  wind  speed 
Uniform  wind  stress  in  x-direction  if  IWIND =0 
Uniform  wind  speed  in  x-direction  if  IWIND = 1 
Uniform  wind  stress  in  y-direction  if  IWIND =0 
Uniform  wind  speed  in  y-direction  if  IWIND  =  1 


DUMMY 
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ISPAC(I),  1=1,10  (1018) 

ISPAC(l)  =  0  ;  Constant  Mannings  n  =  RSPAC(l) 

=  1  ;  Read  Mannings  n  from  File  18 
ISPAC(2  -  3)  =0  ;  Not  used 

ISPAC(4)  =  1  ;  Flag  for  computing  open  boundary  velocities 
ISPAC(5  -  10)  =  0  ;  Not  used 

DUMMY 

JSPAC(I),  1=1,10  (1018) 

JSPAC(l)  =  0  ;  Not  used 

JSPAC(2)  ;  Flag  for  3-D  mode,  quadratic  friction 

=  0  ;  Constant  bottom  friction  factor  =  CTB 

=  1  ;  Bottom  friction  based  on  logarithmic  law 

JSPAC(3)  ;  Flag  for  Coriolis  terms 

=  0  ;  Coriolis  effects  accounted  for 

=  -1  ;  Coriolis  effects  neglected 

JSPAC(4  -  10)  =  0  ;  Not  used 

DUMMY 

RSPAC(I),  1=1,10  (10F8.0) 

RSPAC(l)  ;  Constant  Mannings  n 

RSPAC(2  -  10)  =  0.  ;  Not  used 

DUMMY 

IBTM,  HADD,  HMIN,  HI,  H2,  SSSO,  HMAX  (I8,5F8.0) 

IBTM  ;  Bottom  bathymetry  flag 
IBTM  =  0  ;  Bottom  depth  varies  linearly  from  west  to 
east  of  the  basin 

=  1  ;  Bottom  depth  varies  linearly  from  south  to 

north  of  the  basin 

=  2  ;  Bottom  depth  array  for  cell  center  depths 

read  from  input  file  (FILE  4) 

=  3  ;  Bottom  depth  arrays  HS,  HU,  HV  read  from 

FILE  12 

=  4  ;  Bottom  depths  and  coordinates  of  cell  corners 

read  from  FILE  15  (set  ITRAN=2) 

HADD  ;  A  constant  depth  added  to  the  depth  array 
(cm) 

HMIN  ;  Minimum  water  depth  (cm) 

HI  ;  Bottom  depth  (cm)  along  the  west  or  south 

boundary  of  the  basin  for  IBTM  =  0  or  1 
H2  ;  Bottom  depth  (cm)  along  the  east  or  north 

boundary  of  the  basin  for  IBTM  =  0  or  1 
SSSO  ;  Initial  water  surface  elevation  (cm) 

HMAX  ;  Maximum  water  depth  (cm)  allowed 

DUMMY 

ISMALL,  ISF,  ITB,  ZREFBN,  CTB,  BZl,  ZREFTN,  TZl  (3I8,7F8.0) 
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ISMALL  =  0  ;  Small  amplitude  assumption  is  invoked. 

Surface  elevation  is  not  added  to  the  still 
water  depth  to  obtain  the  total  depth 
=  1  ;  Small  amplitude  assumption  is  not  invoked. 
Surface  elevation  is  added  to  the  still 
water  depth  to  obtain  the  total  depth 
=  0  :  Not  used 


ISF 

ITB 


Bottom  friction  flag 


=  1 
>  1 

ZREFBN 

CTB 

BZl 

ZREFTN 

TZl 


Linear  bottom  friction  for  internal  mode 
Quadratic  bottom  friction  for  internal  mode 
Reference  height  above  bottom  (cm) 
Constant  bottom  drag  coefficient  (typical 
value  0.003) 

;  Bottom  roughness  height  (cm) 

;  Reference  height  at  the  top  (cm) 

;  Constant  surface  roughness  height  (cm) 


DUMMY 

XMAP,  ALXREF,  ALYREF  (10F8.0) 

XMAP  ;  Mapping  factor  that  scales  the  (x,y) 

coordinates  created  by  the  grid  generation 
code  to  the  real  world 

ALXREF  ;  X-reference  length  in  the  computational  plane 
ALYREF  ;  Y-reference  length  in  the  computational  plane 

Note  :  ALXREF  and  ALYREF  are  used  if  ITRAN=  0 

DUMMY 

ITRAN  (1018) 

ITRAN  =  0  ;  Cartesian  grid 

=  1  ;  Curvilinear  grid  created  by  WESCOR.  Cell 
corner  coordinates  rad  from  FILE  15 
=  2  ;  Curvilinear  grid  created  by  WESCORA.  Cell 
corner  coordinates  and  depths  read  from 
FILE15 


DUMMY 

ITBRK(I),  1=1,10  (1018) 

ITBRK(I),  1  =  1,10;  ;  Time  steps  at  which  information  is  written  to 
hot-start  files  (increasing  order) 

DUMMY 

NSTA,  NFREQ,  NSTART  (1018) 

NSTA  ;  Number  of  stations  where  information  is  saved 
for  time  series  plots  of  currents 
NFREQ  ;  Time  step  interval  for  saving  currents 
NSTART  ;  Beginning  time  step  for  saving  currents 

DUMMY 
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IST(K),  JST(K),  STATID(K)  (214, A48) 

If  NSTA  >  0,  NSTA  cards  have  to  be  furnished  below. 

IST(K),  JST(K)  ;  Cell  indices  (I,J)  of  a  station  where 
currents  are  saved 

STATID(K)  ;  Station  description 

DUMMY 

NSTAS,  NFREQS,  NSTRTS  (1018) 

NSTAS  ;  Number  of  stations  where  water  surface 
elevations  are  saved  for  time  series  plots 
NFREQS  ;  Time  step  interval  for  saving  water  surface 
elevations 

NSTRTS  ;  Beginning  time  step  for  saving  water  surface 
elevations 

DUMMY 

ISTS(K),  JSTS(K),  STATS(K)  (214, A48) 

If  NSTAS  >  0,  NSTAS  cards  have  to  be  furnished  below. 

ISTS(K),JSTS(K)  ;  Cell  indices  (I,J)  of  a  station  where  water 
surface  elevations  are  saved 
STATS(K)  ;  Station  description 

DUMMY 

MSTA,  MFREQ,  MSTART  (1018) 

MSTA  ;  Number  of  stations  where  salinity  and 

temperature  information  is  saved  for  time 
series  plots 

MFREQ  ;  Time  step  interval  for  saving  information 
MSTART  ;  Beginning  time  step  for  saving  information 

DUMMY 

ISTSA(K),  JSTSA(K),  STATSA(K)  (2I4,A48) 

If  MSTA  >  0,  MSTA  cards  have  to  be  furnished  below. 
ISTSA(K),JSTSA(K)  ;  Cell  indices  (I,J)  of  a  station  where 
salinity  and  temperature  are  saved 
STATSA(K)  ;  Station  description 

DUMMY 

NRIVER  ;  Number  of  river  boundaries  (2I8,F8. 0,418) 
NRIVER  =0  ;  No  river  boundaries 

<0  ;  River  inflows  are  steady 

>  0  ;  Time  variable  inflows 

If  NRIVER  =  0,  use  the  following  cards 

DUMMY 

DUMMY 

If  NRIVER  >  0,  use  the  following  cards 
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DUMMY 

IJRDIR(K),  IJRROW(K),  IJRSTR(K),  IJREND(K)*  (1018) 


IJRDIR(K)  =  1  ;  River  boundary  is  on  left  (west) 

=  2  ;  River  boundary  is  on  bottom  (south) 

=  3  ;  River  boundary  is  on  right  (east) 

=  4  ;  River  boundary  is  on  top  (north) 

IJRROW(K)  ;  Index  of  the  row  (J)  or  column  (I) 

of  the  river  boundary 

IJRSTR(K)  ;  Starting  I  or  J  index  of  the  river  boundary 
IJREND(K)  ;  Ending  I  or  J  index  of  the  river  boundary 

*  NRIVER  cards  have  to  be  furnished 

If  NRIVER  <  0,  use  the  following  cards 
DUMMY 

IJRDIR(K),  IJRROW(K),  IJRSTR(K),  IJREND(K)‘  (1018) 

IJRDIR(K)  =  1  ;  River  boundary  is  on  left  (west) 

=  2  ;  River  boundary  is  on  bottom  (south) 

=  3  ;  River  boundary  is  on  right  (east) 

=  4  ;  River  boundary  is  on  top  (north) 

IJRROW(K)  ;  Index  of  the  row  (J)  or  column  (I) 

of  the  river  boundary 

IJRSTR(K)  ;  Starting  I  or  J  index  of  the  river  boundary 
IJREND(K)  ;  Ending  I  or  J  index  of  the  river  boundary 

*  I  NRIVER  I  cards  have  to  be  furnished  followed  by  the  cards 
shown  below. 

DUMMY 

ICELL,  JCELL,  QRIVER(K,IJ)*  (2I8,F8.0,4I8) 

ICELL,  JCELL  ;  Coordinates  of  a  cell  (I,J)  where  QRTVER  is 
prescribed 

QRIVER(K,IJ)  ;  Steady  river  inflow 
’Repeat  for  all  the  river  cells,  in  order. 


DUMMY 
NEAR  (1018) 

NEAR  ;  Number  of  interior  thin-wall  barriers 

If  NEAR  =  0,  use  the  following  card 
DUMMY 

If  NEAR  >  0,  use  the  following  cards 
DUMMY 

UEDIR(K),  IJEROW(K),  UESTR(K),  IJEEND(K)*  (1018) 
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IJBDIR(K)  =  1  ;  Barrier  is  in  ^-direction 

=  2  ;  Barrier  is  in  rj-direction 

IJBROW(K)  ;  Index  of  row  (J)  or  column  (I)  of  barrier 

IJBSTR(K)  ;  Starting  I  or  J  index  of  barrier 
IJBEND(K)  ;  Ending  I  or  J  index  of  barrier 
*NBAR  cards  have  to  be  furnished. 


DUMMY 

TIDFNO,  TIDBND  (1018) 

TIDFNO  ;  Number  of  tidal  elevation  tables  entered  as 
input 

TIDBND  ;  Number  of  tidal  elevation  boundaries 
DUMMY 

If  TIDFNO  >  0,  read  the  following  card(s) 

TIDSTR(I),  1=1, TIDFNO  (1018) 

TIDSTR(I)  ;  The  entry  number  in  each  tidal  elevation 
table  corresponding  to  the  starting  time  of 
the  simulation 

DUMMMY 

If  TIDBND  >  0,  TIDBND  cards  of  the  following  format  have  to  be  read. 
IJDIR(I),  UROw’d),  IJSTR(I),  IJEND(I),  TIDTYP(I),  TIDFNl(I), 
TIDFN2(I)  (418, A8, 518) 

IJDIR(I)  =  1  ;  Tidal  boundary  is  on  left  (west) 

=  2  ;  Tidal  boundary  is  on  bottom  (south) 

=  3  ;  Tidal  boundary  is  on  right  (east) 

=  4  ;  Tidal  boundary  is  on  top  (north) 

IJROW(I)  ;  Index  of  the  row  (J)  or  column  (I)  of  the 
tidal  boundary 

IJSTR(I)  ;  Starting  I  or  J  index  of  the  tidal  boundary 
UEND(I)  ;  Ending  I  or  J  index  of  the  tidal  boundary 
TIDTYP(I)  =  "CONSTANT"  ;  Constant  tidal  elevation  between 
IJSTR(I)  and  IJEND(I) 

=  "INTERP  "  ;  Linear  interpolation  of  tidal 
elevation  between  IJSTR(I)  and 
IJEND(I) 

TIDFNl(I)  ;  The  number  of  the  tidal  elevation  table  for 
CONSTANT  or  INTERP  type  of  boundaries 
TIDFN2(I)  ;  The  number  of  the  second  tidal  elevation 

table  used  for  interpolation  on  INTERP  type 
boundaries 


DUMMY 

I,J  (free  format)  ;  Indices  of  a  cell  where  HS  is  reset  to  0. 
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DUMMY 

I,J  (free  format)  ;  Indices  of  a  cell  where  HU  is  reset  to  0. 

DUMMY 

I,J  (free  format)  ;  Indices  of  a  cell  where  HV  is  reset  to  0. 

DUMMY 

I,J,  RDEPTH  (free  format)  ;  Indices  and  depth  (ft)  of  a  cell  where 
HS  is  reset  to  non-zero  value  RDEPTH. 


B10 
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Appendix  C 

List  of  Input  Data  Files 


FILE  13 


River  inflows  are  read  from  FILE  13.  These  data  are  read  first  as  a  time 
line  (DAY  and  HOUR)  formatted  by  218.  Next,  the  (I,J)  location  and  dis¬ 
charge  in  cubic  feet  per  second  for  each  cell  of  each  river  boundary  are  read 
and  formatted  by  (218,  F8.0). 


FILE  14 

Wind  data  are  read  from  FILE  14.  These  data  are  in  the  form  of  time 
(DAY  and  HOUR)  and  the  x  and  y  components  of  the  wind  velocity  in  meters 
per  second  of  each  wind  field  used.  These  data  are  formatted  by 
(2I5,6F10.0). 

FILE  15 


The  (x,y)  coordinates  and  depths  of  the  grid  cell  corners  are  read  from 
FILE  15.  This  file  is  created  from  a  run  of  the  grid  generation  code 
WESCORA  and  a  depth  interpolation  program.  The  first  line  contains  the  file 
name  formatted  as  A80.  The  number  of  comer  points  in  ^  and  r;  are  read 
next  unformatted.  The  coordinates  and  depths  are  read  next  unformatted,  one 
line  per  corner. 

FILE  16 


Tabular  tide  data  are  read  from  FILE  16.  The  first  line  is  the  title  for¬ 
matted  as  A80.  The  tide  data  are  in  the  form  of  time  (MONTH,  DAY, 
YEAR,  HOUR,  MINUTES)  and  the  water  surface  elevations  in  centimeters 
relative  to  selected  datum  for  TIDFNO  points.  These  data  are  formatted  by 
(I2,1X,2I3,1X,2I2,  (T17,8F8.2)). 
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FILE  17 


The  initial  temperature  field  in  degrees  Celsius  is  read  from  FILE  17  by 
format  (10E12.5).  This  file  is  created  from  a  few  observed  values.  The 
resulting  field  is  then  smoothed  in  the  ^  and  rj  directions  several  times  before 
it  is  written  to  FILE  17. 


FILE  18 

A  field  of  Manning’s  n  values  may  be  input  by  format  (20F4.0).  The 
input  values  are  multiplied  by  0.001  in  the  source  code  to  yield  the  actual 
values.  They  are  input  by  rows. 

FILE  19 

Daily  average  equilibrium  temperatures  in  degrees  Celsius  and  surface  heat 
exchange  coefficients  in  units  of  cm/sec  are  read  from  FILE  19.  These  data 
are  in  the  form  of  time  (DAY  and  HOUR),  equilibrium  temperature,  and  heat 
exchange  coefficient.  They  are  formatted  by  (2I5,F10.0,E12.5). 

FILE  74 

The  initial  salinity  field  in  parts  per  thousand  is  read  from  FILE  74  by 
format  (10E12.5).  This  file  is  created  in  the  same  fashion  as  FILE  17. 

FILE  76 


Time-varying  salinity  in  parts  per  thousand  and  temperature  in  degrees 
Celsius  at  tidal  boundaries  are  read  from  FILE  76  if  salinity  and  temperature 
are  to  be  computed.  These  data  are  in  the  form  of  time  (DAY  and  HOUR) 
formatted  by  (215).  Next,  the  (I,J)  location  of  each  tidal  boundary  cell  and 
the  vertical  distribution  of  salinity,  starting  from  the  top  layer  to  the  bottom 
layer  are  read.  These  data  are  followed  by  temperature  data  using  the  same 
format  as  for  the  salinity.  The  format  is  (2I5,11F5.0) 

FILE  78 


Time-varying  temperature  data  at  river  flow  boundaries  are  read  from 
FILE  78  if  temperatures  are  to  be  computed  and  equilibrium  temperatures  are 
not  used  as  river  boundary  temperatures.  These  data  are  in  the  form  of  a  time 
(DAY  and  HOUR)  formatted  by  (215).  Next,  the  (I,J)  location  of  river  flow 
boundary  cells  and  corresponding  temperatures  starting  from  top  layer  to  bot¬ 
tom  layer  are  read.  These  data  are  formatted  by  (2I5,11F6.0). 
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